Analysis for SSFA project: Two factor model

Table of contents

1. Used packages
2. Global settings
3. Load data
4. Model specification
5. Inference
   1. epLsar
   2. R²
   3. Asfc
   4. Smfc
   5. HAsfc9
   6. HAsfc81
6. Summary

Used packages

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import pickle
import arviz as az
import pymc3 as pm

WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

from matplotlib.colors import to_rgb

import scipy.stats as stats 

from IPython.display import display

import matplotlib as mpl

%load_ext autoreload
%autoreload 2

import plotting_lib

Global settings

Output

writeOut = True
outPathPlots = "../plots/statistical_model_two_factors/"
outPathData = "../derived_data/statistical_model_two_factors/"

Plotting

widthMM = 190 
widthInch = widthMM / 25.4
ratio = 0.66666
heigthInch = ratio*widthInch

SMALL_SIZE = 8
MEDIUM_SIZE = 10
BIGGER_SIZE = 12

plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title
sns.set_style("ticks")

dpi = 300

sizes = [SMALL_SIZE,MEDIUM_SIZE,BIGGER_SIZE]

Computing

numSamples = 500
numCores = 10
numTune = 1000
numPredSamples = 2000
random_seed=3651435

Load data

datafile = "../derived_data/preprocessing/preprocessed.dat"

with open(datafile, "rb") as f:
    x1,x2,_,df,dataZ,dictMeanStd,dictTreatment,dictSoftware = pickle.load(f)    

Show that everything is correct:

display(pd.DataFrame.from_dict({'x1':x1,'x2':x2}))

	x1	x2
0	0	5
1	1	5
2	0	5
3	1	5
4	0	5
...	...	...
275	1	9
276	0	9
277	1	9
278	0	9
279	1	9

280 rows × 2 columns

x1 indicates the software used, x2 indicates the treatment applied.

for surfaceParam,(mean,std) in dictMeanStd.items():
    print("Surface parameter {} has mean {} and standard deviation {}".format(surfaceParam,mean,std))

Surface parameter epLsar has mean 0.002866777015225286 and standard deviation 0.0019173233323041528
Surface parameter R² has mean 0.9973197285182144 and standard deviation 0.006745323717352575
Surface parameter Asfc has mean 16.912804592866785 and standard deviation 16.042490777228107
Surface parameter Smfc has mean 2.589874101745179 and standard deviation 10.663178442785044
Surface parameter HAsfc9 has mean 0.31872043020221136 and standard deviation 0.22913790943445264
Surface parameter HAsfc81 has mean 0.5885203613280154 and standard deviation 0.42897734163543366

for k,v in dictTreatment.items():
    print("Number {} encodes treatment {}".format(k,v))

Number 5 encodes treatment Dry bamboo
Number 6 encodes treatment Dry grass
Number 7 encodes treatment Dry lucerne
Number 0 encodes treatment BrushDirt
Number 1 encodes treatment BrushNoDirt
Number 4 encodes treatment Control
Number 10 encodes treatment RubDirt
Number 2 encodes treatment Clover
Number 3 encodes treatment Clover+dust
Number 8 encodes treatment Grass
Number 9 encodes treatment Grass+dust

for k,v in dictSoftware.items():
    print("Number {} encodes software {}".format(k,v))

Number 0 encodes software ConfoMap
Number 1 encodes software Toothfrax

display(dataZ)

	index	TreatmentNumber	SoftwareNumber	DatasetNumber	NameNumber	epLsar_z	R²_z	Asfc_z	Smfc_z	HAsfc9_z	HAsfc81_z
0	0	5	0	0	116	0.896441	0.202280	-0.548753	-0.173680	-0.872185	-0.512580
1	1	5	1	0	116	0.967089	0.332122	-0.410913	-0.231700	-0.799735	-0.528438
2	2	5	0	0	117	1.663582	0.205805	-0.406276	-0.184263	-0.625739	-0.636413
3	3	5	1	0	117	1.559060	0.318335	-0.231484	-0.231700	-0.652395	-0.762985
4	4	5	0	0	118	1.235447	0.140997	-0.524577	-0.187419	-0.425552	-0.460441
...	...	...	...	...	...	...	...	...	...	...	...
275	275	9	1	2	51	0.812186	0.302472	-0.940814	-0.110773	2.259947	1.219611
276	276	9	0	2	52	0.279516	0.134624	-0.883003	-0.199355	1.616592	2.457881
277	277	9	1	2	52	0.141981	0.358659	-0.882314	-0.230373	2.779893	2.898056
278	278	9	0	2	53	-0.859089	0.269132	-0.964638	-0.118959	0.316720	0.540611
279	279	9	1	2	53	-1.128540	0.376153	-0.964978	-0.204577	0.119539	0.586105

280 rows × 11 columns

display(df)

	Dataset	Name	Software	Diet	Treatment	Before.after	epLsar	R²	Asfc	Smfc	HAsfc9	HAsfc81	NewEplsar	TreatmentNumber	SoftwareNumber	DatasetNumber	NameNumber
0	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_1	ConfoMap	Dry bamboo	Dry bamboo	NaN	0.004586	0.998684	8.109445	0.737898	0.118870	0.368635	0.019529	5	0	0	116
1	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_1	Toothfrax	Dry bamboo	Dry bamboo	NaN	0.004721	0.999560	10.320730	0.119219	0.135471	0.361833	NaN	5	1	0	116
2	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_2	ConfoMap	Dry bamboo	Dry bamboo	NaN	0.006056	0.998708	10.395128	0.625040	0.175340	0.315513	0.020162	5	0	0	117
3	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_2	Toothfrax	Dry bamboo	Dry bamboo	NaN	0.005856	0.999467	13.199232	0.119219	0.169232	0.261217	NaN	5	1	0	117
4	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_3	ConfoMap	Dry bamboo	Dry bamboo	NaN	0.005236	0.998271	8.497286	0.591396	0.221210	0.391002	0.019804	5	0	0	118
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
275	Sheeps	L8-Ovis-90730-lm2sin-a	Toothfrax	Grass+dust	Grass+dust	NaN	0.004424	0.999360	1.819802	1.408678	0.836560	1.111706	NaN	9	1	2	51
276	Sheeps	L8-Ovis-90764-lm2sin-a	ConfoMap	Grass+dust	Grass+dust	NaN	0.003403	0.998228	2.747241	0.464115	0.689143	1.642896	0.018978	9	0	2	52
277	Sheeps	L8-Ovis-90764-lm2sin-a	Toothfrax	Grass+dust	Grass+dust	NaN	0.003139	0.999739	2.758297	0.133366	0.955699	1.831721	NaN	9	1	2	52
278	Sheeps	L8-Ovis-90814-lm2sin-a	ConfoMap	Grass+dust	Grass+dust	NaN	0.001220	0.999135	1.437607	1.321398	0.391293	0.820430	0.017498	9	0	2	53
279	Sheeps	L8-Ovis-90814-lm2sin-a	Toothfrax	Grass+dust	Grass+dust	NaN	0.000703	0.999857	1.432148	0.408433	0.346111	0.839946	NaN	9	1	2	53

280 rows × 17 columns

Model specification

class TwoFactorModel(pm.Model):
    
    """
    Compute params of priors and hyperpriors.
    """
    def getParams(self,x1,x2,y):
        # get lengths
        Nx1Lvl = np.unique(x1).size
        Nx2Lvl = np.unique(x2).size
                        
        dims = (Nx1Lvl, Nx2Lvl)
        
        ### get standard deviations
        
        # convert to pandas dataframe to use their logic
        df = pd.DataFrame.from_dict({'x1':x1,'x2':x2,'y':y})
        
        s1 = df.groupby('x1').std()['y'].max()
        s2 = df.groupby('x2').std()['y'].max()
                
        stdSingle = (s1, s2)
        
        prefac = 0.05
        s12 = prefac * np.linalg.norm([s1,s2])
        
        stdMulti = (s12)
        
        return (dims, stdSingle, stdMulti)
    
    def printParams(self,x1,x2,y):
        dims, stdSingle, stdMulti = self.getParams(x1,x2,y)
        Nx1Lvl, Nx2Lvl = dims
        s1, s2 = stdSingle
        s12 = stdMulti
        
        print("The number of levels of the x variables are {}".format(dims))
        print("The standard deviations used for the beta priors are  {}".format(stdSingle))
        print("The standard deviations used for the M12 priors are {}".format(stdMulti))
    
    def __init__(self,name,x1,x2,y,model=None):
        
        # call super's init first, passing model and name
        super().__init__(name, model)
        
        # get parameter of hyperpriors
        dims, stdSingle, stdMulti = self.getParams(x1,x2,y)
        Nx1Lvl, Nx2Lvl = dims
        s1, s2 = stdSingle
        s12 = stdMulti
        
        ### hyperpriors ### 
        # observation hyperpriors
        lamY = 1/30.
        muGamma = 0.5
        sigmaGamma = 2.
        
        # prediction hyperpriors
        sigma0 = pm.HalfNormal('sigma0',sd=1)
        sigma1 = pm.HalfNormal('sigma1',sd=s1, shape=Nx1Lvl)       
        sigma2 = pm.HalfNormal('sigma2',sd=s2, shape=Nx2Lvl)
        beta2 = (np.sqrt(6)*sigma2)/(np.pi)
        
        mu_b0 = pm.Normal('mu_b0', mu=0., sd=1)
        mu_b1 = pm.Normal('mu_b1', mu=0., sd=1, shape=Nx1Lvl)        
        mu_b2 = pm.Normal('mu_b2', mu=0., sd=1, shape=Nx2Lvl)       
        
        sigma12 = pm.HalfNormal('sigma12',sd=s12)
                               
        ### priors ### 
        # observation priors        
        nuY = pm.Exponential('nuY',lam=lamY)
        sigmaY = pm.Gamma('sigmaY',mu=muGamma, sigma=sigmaGamma)
        
        # prediction priors
        b0_dist = pm.Normal('b0_dist', mu=0, sd=1)
        b0 = pm.Deterministic("b0", mu_b0 + b0_dist * sigma0)
       
        b1_dist = pm.Normal('b1_dist', mu=0, sd=1)
        b1 = pm.Deterministic("b1", mu_b1 + b1_dist * sigma1)
                        
        b2_beta = pm.HalfNormal('b2_beta', sd=beta2, shape=Nx2Lvl)
        b2_dist = pm.Gumbel('b2_dist', mu=0, beta=1)
        b2 = pm.Deterministic("b2", mu_b2 + b2_beta * b2_dist)
        
        mu_M12 = pm.Normal('mu_M12', mu=0., sd=1, shape=[Nx1Lvl, Nx2Lvl])
        M12_dist = pm.Normal('M12_dist', mu=0, sd=1)
        M12 = pm.Deterministic("M12", mu_M12 + M12_dist * sigma12)
        
        #### prediction ###         
        mu = pm.Deterministic('mu',b0 + b1[x1]+ b2[x2] +  M12[x1,x2] )
                                        
        ### observation ### 
        y = pm.StudentT('y',nu = nuY, mu=mu, sd=sigmaY, observed=y)

Inference

epLsar

with pm.Model() as model:
    epLsarModel = TwoFactorModel('epLsar',x1,x2,dataZ.epLsar_z.values)

Verify model settings

epLsarModel.printParams(x1,x2,dataZ.epLsar_z.values)

The number of levels of the x variables are (2, 11)
The standard deviations used for the beta priors are  (1.0100483277420549, 1.278487256789849)
The standard deviations used for the M12 priors are 0.08146666941376325

try:
    graph_epLsar = pm.model_to_graphviz(epLsarModel)    
except:
    graph_epLsar = "Could not make graph"
graph_epLsar

Check prior choice

with epLsarModel as model:
    prior_pred_epLsar = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_epLsar,dataZ.epLsar_z.values,'epLsar')

Prior choice is as intended: Broad over the data range.

Sampling

with epLsarModel as model:
    trace_epLsar = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.95,random_seed=random_seed)
    #fit_epLsar = pm.fit(random_seed=random_seed)
    #trace_epLsar = fit_epLsar.sample(draws=numSamples)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [epLsar_M12_dist, epLsar_mu_M12, epLsar_b2_dist, epLsar_b2_beta, epLsar_b1_dist, epLsar_b0_dist, epLsar_sigmaY, epLsar_nuY, epLsar_sigma12, epLsar_mu_b2, epLsar_mu_b1, epLsar_mu_b0, epLsar_sigma2, epLsar_sigma1, epLsar_sigma0]

100.00% [15000/15000 03:01<00:00 Sampling 10 chains, 0 divergences]

Sampling 10 chains for 1_000 tune and 500 draw iterations (10_000 + 5_000 draws total) took 183 seconds.

with epLsarModel as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('epLsar'), 'wb') as buff:
            pickle.dump({'model':epLsarModel, 'trace': trace_epLsar}, buff)            

Save for later comparison

if writeOut:
    np.save('../derived_data/statistical_model_two_factors/epLsar_oldb1', trace_epLsar['epLsar_b1'])
    np.save('../derived_data/statistical_model_two_factors/epLsar_oldb2', trace_epLsar['epLsar_b2'])
    np.save('../derived_data/statistical_model_two_factors/epLsar_oldM12', trace_epLsar['epLsar_M12'])

Check sampling

with epLsarModel as model:
    dataTrace_epLsar = az.from_pymc3(trace=trace_epLsar)

pm.summary(dataTrace_epLsar,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
epLsar_mu_b0	-0.11	0.81	-1.67	1.48	0.01	0.01	5029.0	2319.0	5029.0	3627.0	1.0
epLsar_mu_b1[0]	0.00	0.73	-1.48	1.36	0.01	0.01	4766.0	3057.0	4763.0	3988.0	1.0
epLsar_mu_b1[1]	-0.07	0.71	-1.47	1.34	0.01	0.01	4575.0	2861.0	4587.0	3606.0	1.0
epLsar_mu_b2[0]	-0.25	0.66	-1.63	0.94	0.01	0.01	3705.0	2967.0	3733.0	4054.0	1.0
epLsar_mu_b2[1]	-0.36	0.67	-1.64	0.98	0.01	0.01	4093.0	3216.0	4098.0	3736.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
epLsar_mu[275]	0.66	0.27	0.12	1.16	0.00	0.00	5124.0	4919.0	5132.0	4303.0	1.0
epLsar_mu[276]	0.68	0.31	0.09	1.30	0.00	0.00	5093.0	5024.0	5103.0	3850.0	1.0
epLsar_mu[277]	0.66	0.27	0.12	1.16	0.00	0.00	5124.0	4919.0	5132.0	4303.0	1.0
epLsar_mu[278]	0.68	0.31	0.09	1.30	0.00	0.00	5093.0	5024.0	5103.0	3850.0	1.0
epLsar_mu[279]	0.66	0.27	0.12	1.16	0.00	0.00	5124.0	4919.0	5132.0	4303.0	1.0

384 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_epLsar,dataTrace_epLsar,'epLsar')

/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  warnings.warn(

with epLsarModel as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_epLsar,'epLsar')

with epLsarModel as model:
    plotting_lib.pm.energyplot(trace_epLsar)

Posterior predictive distribution

with epLsarModel as model:
    posterior_pred_epLsar = pm.sample_posterior_predictive(trace_epLsar,samples=numPredSamples,random_seed=random_seed)

/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  warnings.warn(

100.00% [2000/2000 00:02<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_epLsar,posterior_pred_epLsar,dataZ.epLsar_z.values,'epLsar')

Level plots

plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_epLsar,'epLsar',x1,x2)

plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_epLsar,'epLsar',x1,x2)

Posterior and contrasts

df_hdi_epLsar = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_epLsar,'epLsar',x1,x2)

df_hdi_epLsar

	Treatment_i	Treatment_j	hdi_ConfoMap_2.5%	hdi_ConfoMap_97.5%	isSignificant_on_ConfoMap	hdi_Toothfrax_2.5%	hdi_Toothfrax_97.5%	isSignificant_on_Toothfrax
0	Dry grass	Dry bamboo	-0.004200	-0.002884	True	-0.003967	-0.002636	True
1	Dry lucerne	Dry bamboo	-0.003305	-0.001994	True	-0.003411	-0.002161	True
2	Dry lucerne	Dry grass	0.000228	0.001524	True	-0.000167	0.001116	False
3	BrushNoDirt	BrushDirt	-0.001227	0.001107	False	-0.001528	0.000527	False
4	Control	BrushDirt	-0.002046	0.000359	False	-0.001954	0.000235	False
5	Control	BrushNoDirt	-0.001872	0.000384	False	-0.001414	0.000676	False
6	RubDirt	BrushDirt	-0.001636	0.000627	False	-0.001497	0.000683	False
7	RubDirt	BrushNoDirt	-0.001442	0.000661	False	-0.000890	0.001130	False
8	RubDirt	Control	-0.000753	0.001508	False	-0.000598	0.001547	False
9	Clover+dust	Clover	-0.001140	0.001047	False	-0.000746	0.001576	False
10	Grass	Clover	-0.000153	0.002973	False	0.000293	0.003795	True
11	Grass	Clover+dust	-0.000151	0.002972	False	-0.000175	0.003389	False
12	Grass+dust	Clover	0.000533	0.003296	True	0.000819	0.003383	True
13	Grass+dust	Clover+dust	0.000445	0.003307	True	0.000379	0.002978	True
14	Grass+dust	Grass	-0.001303	0.002230	False	-0.001721	0.001958	False

if writeOut:
    df_hdi_epLsar.to_csv(outPathData+ 'hdi_{}.csv'.format('epLsar'))

R²

with pm.Model() as model:
    RsquaredModel = TwoFactorModel('R²',x1,x2,dataZ["R²_z"].values)

Verify model settings

RsquaredModel.printParams(x1,x2,dataZ["R²_z"].values)

The number of levels of the x variables are (2, 11)
The standard deviations used for the beta priors are  (1.3433087377122408, 2.3171005331515255)
The standard deviations used for the M12 priors are 0.13391632877981252

pm.model_to_graphviz(RsquaredModel)

Check prior choice

with RsquaredModel as model:
    prior_pred_Rsquared = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Rsquared,dataZ["R²_z"].values,'R²')

Sampling

with RsquaredModel as model:
    trace_Rsquared = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [R²_M12_dist, R²_mu_M12, R²_b2_dist, R²_b2_beta, R²_b1_dist, R²_b0_dist, R²_sigmaY, R²_nuY, R²_sigma12, R²_mu_b2, R²_mu_b1, R²_mu_b0, R²_sigma2, R²_sigma1, R²_sigma0]

100.00% [15000/15000 1:06:10<00:00 Sampling 10 chains, 0 divergences]

Sampling 10 chains for 1_000 tune and 500 draw iterations (10_000 + 5_000 draws total) took 3971 seconds.

with RsquaredModel as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('R²'), 'wb') as buff:
            pickle.dump({'model': RsquaredModel, 'trace': trace_Rsquared}, buff)

Check sampling

with RsquaredModel as model:
    dataTrace_Rsquared = az.from_pymc3(trace=trace_Rsquared)

pm.summary(dataTrace_Rsquared,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
R²_mu_b0	0.07	0.83	-1.55	1.64	0.01	0.01	7950.0	2327.0	7929.0	4058.0	1.0
R²_mu_b1[0]	-0.05	0.73	-1.49	1.35	0.01	0.01	6903.0	2379.0	6902.0	3964.0	1.0
R²_mu_b1[1]	0.11	0.72	-1.23	1.58	0.01	0.01	7161.0	2924.0	7164.0	4300.0	1.0
R²_mu_b2[0]	-0.01	0.66	-1.24	1.29	0.01	0.01	5941.0	2544.0	5938.0	3984.0	1.0
R²_mu_b2[1]	-0.05	0.66	-1.31	1.29	0.01	0.01	6122.0	2330.0	6086.0	3681.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
R²_mu[275]	0.34	0.01	0.31	0.36	0.00	0.00	4973.0	4968.0	4975.0	4789.0	1.0
R²_mu[276]	0.11	0.10	-0.03	0.26	0.00	0.00	4537.0	4387.0	4704.0	4535.0	1.0
R²_mu[277]	0.34	0.01	0.31	0.36	0.00	0.00	4973.0	4968.0	4975.0	4789.0	1.0
R²_mu[278]	0.11	0.10	-0.03	0.26	0.00	0.00	4537.0	4387.0	4704.0	4535.0	1.0
R²_mu[279]	0.34	0.01	0.31	0.36	0.00	0.00	4973.0	4968.0	4975.0	4789.0	1.0

384 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Rsquared,dataTrace_Rsquared,'R²')

/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  warnings.warn(

with RsquaredModel as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Rsquared,'R²')

with RsquaredModel as model:
    plotting_lib.pm.energyplot(trace_Rsquared)

Posterior predictive distribution

with RsquaredModel as model:
    posterior_pred_Rsquared = pm.sample_posterior_predictive(trace_Rsquared,samples=numPredSamples,random_seed=random_seed)

/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  warnings.warn(

100.00% [2000/2000 00:02<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Rsquared,posterior_pred_Rsquared,dataZ["R²_z"].values,'R²')

Compare prior and posterior for model parameters

with RsquaredModel as model:
    pm_data_Rsquared = az.from_pymc3(trace=trace_Rsquared,prior=prior_pred_Rsquared,posterior_predictive=posterior_pred_Rsquared)

arviz.data.io_pymc3 - WARNING - posterior predictive variable R²_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_Rsquared,'R²')

plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Rsquared,'R²',x1,x2)

plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Rsquared,'R²',x1,x2)

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_Rsquared,'R²')

df_hdi_R = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Rsquared,'R²',x1,x2)

/home/bob/.local/lib/python3.8/site-packages/seaborn/distributions.py:305: UserWarning: Dataset has 0 variance; skipping density estimate.
  warnings.warn(msg, UserWarning)

df_hdi_R

	Treatment_i	Treatment_j	hdi_ConfoMap_2.5%	hdi_ConfoMap_97.5%	isSignificant_on_ConfoMap	hdi_Toothfrax_2.5%	hdi_Toothfrax_97.5%	isSignificant_on_Toothfrax
0	Dry grass	Dry bamboo	-0.000285	0.000175	False	-0.000257	0.000188	False
1	Dry lucerne	Dry bamboo	-0.000116	0.000248	False	-0.000326	0.000121	False
2	Dry lucerne	Dry grass	-0.000115	0.000330	False	-0.000349	0.000179	False
3	BrushNoDirt	BrushDirt	-0.002096	-0.000672	True	-0.000420	0.000635	False
4	Control	BrushDirt	-0.001249	0.000041	False	-0.000416	0.000408	False
5	Control	BrushNoDirt	-0.000080	0.001656	False	-0.000801	0.000387	False
6	RubDirt	BrushDirt	-0.001154	0.000799	False	-0.000185	0.000635	False
7	RubDirt	BrushNoDirt	0.000181	0.002386	True	-0.000556	0.000619	False
8	RubDirt	Control	-0.000680	0.001489	False	-0.000252	0.000722	False
9	Clover+dust	Clover	-0.001005	-0.000290	True	-0.000138	0.000401	False
10	Grass	Clover	-0.001250	-0.000102	True	-0.000719	-0.000045	True
11	Grass	Clover+dust	-0.000639	0.000580	False	-0.000792	-0.000199	True
12	Grass+dust	Clover	-0.001970	0.000084	False	-0.000435	0.000140	False
13	Grass+dust	Clover+dust	-0.001314	0.000824	False	-0.000515	-0.000045	True
14	Grass+dust	Grass	-0.001536	0.000874	False	-0.000091	0.000540	False

if writeOut:
    df_hdi_R.to_csv(outPathData+ 'hdi_{}.csv'.format('R²'))

Asfc

with pm.Model() as model:
    AsfcModel = TwoFactorModel('Asfc',x1,x2,dataZ["Asfc_z"].values)

Verify model settings

AsfcModel.printParams(x1,x2,dataZ["Asfc_z"].values)

The number of levels of the x variables are (2, 11)
The standard deviations used for the beta priors are  (1.0620939450667206, 0.7885965314025212)
The standard deviations used for the M12 priors are 0.06614242279897746

pm.model_to_graphviz(AsfcModel)

Check prior choice

with AsfcModel as model:
    prior_pred_Asfc = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Asfc,dataZ["Asfc_z"].values,'Asfc')

Prior choice is as intended: Broad over the data range.

Sampling

with AsfcModel as model:
    trace_Asfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [Asfc_M12_dist, Asfc_mu_M12, Asfc_b2_dist, Asfc_b2_beta, Asfc_b1_dist, Asfc_b0_dist, Asfc_sigmaY, Asfc_nuY, Asfc_sigma12, Asfc_mu_b2, Asfc_mu_b1, Asfc_mu_b0, Asfc_sigma2, Asfc_sigma1, Asfc_sigma0]

100.00% [15000/15000 17:10<00:00 Sampling 10 chains, 0 divergences]

Sampling 10 chains for 1_000 tune and 500 draw iterations (10_000 + 5_000 draws total) took 1032 seconds.

with AsfcModel as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('Asfc'), 'wb') as buff:
            pickle.dump({'model': AsfcModel, 'trace': trace_Asfc}, buff)

Check sampling

with AsfcModel as model:
    dataTrace_Asfc = az.from_pymc3(trace=trace_Asfc)

pm.summary(dataTrace_Asfc,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
Asfc_mu_b0	0.03	0.79	-1.54	1.57	0.01	0.01	7171.0	2541.0	7170.0	3979.0	1.0
Asfc_mu_b1[0]	-0.08	0.70	-1.44	1.28	0.01	0.01	6377.0	2385.0	6384.0	3933.0	1.0
Asfc_mu_b1[1]	0.11	0.71	-1.33	1.47	0.01	0.01	6789.0	2647.0	6777.0	4089.0	1.0
Asfc_mu_b2[0]	0.67	0.68	-0.75	1.92	0.01	0.01	6024.0	4059.0	6090.0	3494.0	1.0
Asfc_mu_b2[1]	0.80	0.68	-0.50	2.16	0.01	0.01	6900.0	4532.0	6872.0	3974.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
Asfc_mu[275]	-0.94	0.05	-1.04	-0.84	0.00	0.00	4903.0	4900.0	4919.0	4492.0	1.0
Asfc_mu[276]	-0.94	0.05	-1.04	-0.84	0.00	0.00	5312.0	5312.0	5324.0	4413.0	1.0
Asfc_mu[277]	-0.94	0.05	-1.04	-0.84	0.00	0.00	4903.0	4900.0	4919.0	4492.0	1.0
Asfc_mu[278]	-0.94	0.05	-1.04	-0.84	0.00	0.00	5312.0	5312.0	5324.0	4413.0	1.0
Asfc_mu[279]	-0.94	0.05	-1.04	-0.84	0.00	0.00	4903.0	4900.0	4919.0	4492.0	1.0

384 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Asfc,dataTrace_Asfc,'Asfc')

/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  warnings.warn(

with AsfcModel as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Asfc,'Asfc')

with AsfcModel as model:
    plotting_lib.pm.energyplot(trace_Asfc)

Posterior predictive distribution

with AsfcModel as model:
    posterior_pred_Asfc = pm.sample_posterior_predictive(trace_Asfc,samples=numPredSamples,random_seed=random_seed)

/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  warnings.warn(

100.00% [2000/2000 00:03<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Asfc,posterior_pred_Asfc,dataZ["Asfc_z"].values,'Asfc')

plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Asfc,'Asfc',x1,x2)

plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Asfc,'Asfc',x1,x2)

Compare prior and posterior for model parameters

with AsfcModel as model:
    pm_data_Asfc = az.from_pymc3(trace=trace_Asfc,prior=prior_pred_Asfc,posterior_predictive=posterior_pred_Asfc)

arviz.data.io_pymc3 - WARNING - posterior predictive variable Asfc_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_Asfc,'Asfc')

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_Asfc,'Asfc')

df_hdi_Asfc = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Asfc,'Asfc',x1,x2)

df_hdi_Asfc

	Treatment_i	Treatment_j	hdi_ConfoMap_2.5%	hdi_ConfoMap_97.5%	isSignificant_on_ConfoMap	hdi_Toothfrax_2.5%	hdi_Toothfrax_97.5%	isSignificant_on_Toothfrax
0	Dry grass	Dry bamboo	-3.939803	0.935695	False	-4.447570	0.965258	False
1	Dry lucerne	Dry bamboo	-4.031412	0.625569	False	-5.003380	0.239099	False
2	Dry lucerne	Dry grass	-2.410119	1.424301	False	-2.515884	1.276464	False
3	BrushNoDirt	BrushDirt	0.013059	7.915624	True	-2.631889	6.887819	False
4	Control	BrushDirt	3.910899	18.434847	True	1.604132	14.612596	True
5	Control	BrushNoDirt	0.022112	14.415237	True	-0.616464	13.182578	False
6	RubDirt	BrushDirt	-0.259987	20.831039	False	1.502426	13.522475	True
7	RubDirt	BrushNoDirt	-4.336062	16.137641	False	-1.231718	11.842114	False
8	RubDirt	Control	-14.247123	11.427003	False	-8.595219	7.100929	False
9	Clover+dust	Clover	-4.902063	2.035175	False	-4.562186	1.844004	False
10	Grass	Clover	-5.948810	1.110713	False	-5.876442	0.688654	False
11	Grass	Clover+dust	-3.539383	1.255674	False	-3.501566	1.339794	False
12	Grass+dust	Clover	-6.781348	0.056399	False	-6.135492	0.254948	False
13	Grass+dust	Clover+dust	-4.430594	0.335530	False	-3.919102	0.673998	False
14	Grass+dust	Grass	-3.353905	1.617864	False	-3.044589	1.686165	False

if writeOut:
    df_hdi_Asfc.to_csv(outPathData+ 'hdi_{}.csv'.format('Asfc'))

Smfc

with pm.Model() as model:
    SmfcModel = TwoFactorModel('Smfc',x1,x2,dataZ.Smfc_z.values)

Verify model settings

SmfcModel.printParams(x1,x2,dataZ.Smfc_z.values)

The number of levels of the x variables are (2, 11)
The standard deviations used for the beta priors are  (1.190539275484547, 2.890262579386322)
The standard deviations used for the M12 priors are 0.15629300643560595

pm.model_to_graphviz(SmfcModel)

Check prior choice

with SmfcModel as model:
    prior_pred_Smfc = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Smfc,dataZ.Smfc_z.values,'Smfc')

Prior choice is as intended: Broad over the data range.

Sampling

with SmfcModel as model:
    trace_Smfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [Smfc_M12_dist, Smfc_mu_M12, Smfc_b2_dist, Smfc_b2_beta, Smfc_b1_dist, Smfc_b0_dist, Smfc_sigmaY, Smfc_nuY, Smfc_sigma12, Smfc_mu_b2, Smfc_mu_b1, Smfc_mu_b0, Smfc_sigma2, Smfc_sigma1, Smfc_sigma0]

15.05% [2257/15000 18:47<1:46:07 Sampling 10 chains, 0 divergences]

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
~/.local/lib/python3.8/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs)
   1485             with sampler:
-> 1486                 for draw in sampler:
   1487                     trace = traces[draw.chain - chain]

~/.local/lib/python3.8/site-packages/pymc3/parallel_sampling.py in __iter__(self)
    491         while self._active:
--> 492             draw = ProcessAdapter.recv_draw(self._active)
    493             proc, is_last, draw, tuning, stats, warns = draw

~/.local/lib/python3.8/site-packages/pymc3/parallel_sampling.py in recv_draw(processes, timeout)
    351         pipes = [proc._msg_pipe for proc in processes]
--> 352         ready = multiprocessing.connection.wait(pipes)
    353         if not ready:

/usr/lib/python3.8/multiprocessing/connection.py in wait(object_list, timeout)
    930             while True:
--> 931                 ready = selector.select(timeout)
    932                 if ready:

/usr/lib/python3.8/selectors.py in select(self, timeout)
    414         try:
--> 415             fd_event_list = self._selector.poll(timeout)
    416         except InterruptedError:

KeyboardInterrupt: 

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-90-4de91012caf3> in <module>
      1 with SmfcModel as model:
----> 2     trace_Smfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

~/.local/lib/python3.8/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs)
    543         _print_step_hierarchy(step)
    544         try:
--> 545             trace = _mp_sample(**sample_args, **parallel_args)
    546         except pickle.PickleError:
    547             _log.warning("Could not pickle model, sampling singlethreaded.")

~/.local/lib/python3.8/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs)
   1510     except KeyboardInterrupt:
   1511         if discard_tuned_samples:
-> 1512             traces, length = _choose_chains(traces, tune)
   1513         else:
   1514             traces, length = _choose_chains(traces, 0)

~/.local/lib/python3.8/site-packages/pymc3/sampling.py in _choose_chains(traces, tune)
   1528     lengths = [max(0, len(trace) - tune) for trace in traces]
   1529     if not sum(lengths):
-> 1530         raise ValueError("Not enough samples to build a trace.")
   1531 
   1532     idxs = np.argsort(lengths)[::-1]

ValueError: Not enough samples to build a trace.

Analysis stopped here because sampling did not converge. As the plot shows, some data points are very far away from the others, which would require the analysis to be based on more heavy-tailed distributions.

HAsfc9

with pm.Model() as model:
    HAsfc9Model = TwoFactorModel('HAsfc9',x1,x2,dataZ["HAsfc9_z"].values)

Verify model settings

HAsfc9Model.printParams(x1,x2,dataZ["HAsfc9_z"].values)

The number of levels of the x variables are (2, 11)
The standard deviations used for the beta priors are  (1.0540496136136044, 2.005676747692769)
The standard deviations used for the M12 priors are 0.11328900878057889

pm.model_to_graphviz(HAsfc9Model)

Check prior choice

with HAsfc9Model as model:
    prior_pred_HAsfc9 = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc9,dataZ["HAsfc9_z"].values,'HAsfc9')

Prior choice is as intended: Broad over the data range.

Sampling

with HAsfc9Model as model:
    trace_HAsfc9 = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [HAsfc9_M12_dist, HAsfc9_mu_M12, HAsfc9_b2_dist, HAsfc9_b2_beta, HAsfc9_b1_dist, HAsfc9_b0_dist, HAsfc9_sigmaY, HAsfc9_nuY, HAsfc9_sigma12, HAsfc9_mu_b2, HAsfc9_mu_b1, HAsfc9_mu_b0, HAsfc9_sigma2, HAsfc9_sigma1, HAsfc9_sigma0]

100.00% [15000/15000 06:51<00:00 Sampling 10 chains, 0 divergences]

Sampling 10 chains for 1_000 tune and 500 draw iterations (10_000 + 5_000 draws total) took 412 seconds.

with HAsfc9Model as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('HAsfc9'), 'wb') as buff:
            pickle.dump({'model': HAsfc9Model, 'trace': trace_HAsfc9}, buff)

Check sampling

with HAsfc9Model as model:
    dataTrace_HAsfc9 = az.from_pymc3(trace=trace_HAsfc9)

pm.summary(dataTrace_HAsfc9,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
HAsfc9_mu_b0	-0.07	0.80	-1.65	1.43	0.01	0.01	5695.0	2359.0	5673.0	4066.0	1.0
HAsfc9_mu_b1[0]	0.04	0.72	-1.34	1.51	0.01	0.01	5050.0	2946.0	5057.0	3783.0	1.0
HAsfc9_mu_b1[1]	-0.10	0.72	-1.49	1.35	0.01	0.01	4834.0	2979.0	4834.0	3973.0	1.0
HAsfc9_mu_b2[0]	-0.26	0.66	-1.56	1.03	0.01	0.01	4119.0	2839.0	4127.0	3512.0	1.0
HAsfc9_mu_b2[1]	-0.20	0.67	-1.57	1.08	0.01	0.01	3738.0	3165.0	3749.0	3996.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
HAsfc9_mu[275]	-0.03	0.19	-0.41	0.34	0.00	0.00	4547.0	4013.0	4560.0	4624.0	1.0
HAsfc9_mu[276]	0.06	0.24	-0.39	0.52	0.00	0.00	4766.0	2887.0	4774.0	4037.0	1.0
HAsfc9_mu[277]	-0.03	0.19	-0.41	0.34	0.00	0.00	4547.0	4013.0	4560.0	4624.0	1.0
HAsfc9_mu[278]	0.06	0.24	-0.39	0.52	0.00	0.00	4766.0	2887.0	4774.0	4037.0	1.0
HAsfc9_mu[279]	-0.03	0.19	-0.41	0.34	0.00	0.00	4547.0	4013.0	4560.0	4624.0	1.0

384 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc9,dataTrace_HAsfc9,'HAsfc9')

/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  warnings.warn(

with HAsfc9Model as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc9,'HAsfc9')

with HAsfc9Model as model:
    plotting_lib.pm.energyplot(trace_HAsfc9)

Posterior predictive distribution

with HAsfc9Model as model:
    posterior_pred_HAsfc9 = pm.sample_posterior_predictive(trace_HAsfc9,samples=numPredSamples,random_seed=random_seed)

/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  warnings.warn(

100.00% [2000/2000 00:02<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc9,posterior_pred_HAsfc9,dataZ["HAsfc9_z"].values,'HAsfc9')

plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc9,'HAsfc9',x1,x2)

plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc9,'HAsfc9',x1,x2)

Compare prior and posterior for model parameters

with HAsfc9Model as model:
    pm_data_HAsfc9 = az.from_pymc3(trace=trace_HAsfc9,prior=prior_pred_HAsfc9,posterior_predictive=posterior_pred_HAsfc9)

arviz.data.io_pymc3 - WARNING - posterior predictive variable HAsfc9_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc9,'HAsfc9')

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc9,'HAsfc9')

df_hdi_HAsfc9 = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc9,'HAsfc9',x1,x2)

df_hdi_HAsfc9

	Treatment_i	Treatment_j	hdi_ConfoMap_2.5%	hdi_ConfoMap_97.5%	isSignificant_on_ConfoMap	hdi_Toothfrax_2.5%	hdi_Toothfrax_97.5%	isSignificant_on_Toothfrax
0	Dry grass	Dry bamboo	0.000670	0.138311	True	0.002941	0.143668	True
1	Dry lucerne	Dry bamboo	-0.035688	0.084977	False	-0.042002	0.076006	False
2	Dry lucerne	Dry grass	-0.117764	0.031872	False	-0.125602	0.021941	False
3	BrushNoDirt	BrushDirt	-0.231361	0.146793	False	-0.063532	0.221075	False
4	Control	BrushDirt	-0.159309	0.177833	False	-0.020488	0.239702	False
5	Control	BrushNoDirt	-0.118868	0.215704	False	-0.121261	0.192874	False
6	RubDirt	BrushDirt	-0.223620	0.119671	False	-0.093554	0.146212	False
7	RubDirt	BrushNoDirt	-0.184998	0.151141	False	-0.204337	0.093804	False
8	RubDirt	Control	-0.211662	0.080935	False	-0.212686	0.056618	False
9	Clover+dust	Clover	-0.147903	0.136656	False	-0.081257	0.155239	False
10	Grass	Clover	-0.251859	0.071005	False	-0.126162	0.116917	False
11	Grass	Clover+dust	-0.231351	0.051562	False	-0.160442	0.079536	False
12	Grass+dust	Clover	-0.310090	0.004148	False	-0.202760	0.043221	False
13	Grass+dust	Clover+dust	-0.272959	-0.010159	True	-0.240244	-0.005423	True
14	Grass+dust	Grass	-0.211711	0.091591	False	-0.206979	0.039560	False

if writeOut:
    df_hdi_HAsfc9.to_csv(outPathData+ 'hdi_{}.csv'.format('HAsfc9'))

HAsfc81

with pm.Model() as model:
    HAsfc81Model = TwoFactorModel('HAsfc81',x1,x2,dataZ["HAsfc81_z"].values)

Verify model settings

HAsfc81Model.printParams(x1,x2,dataZ["HAsfc81_z"].values)

The number of levels of the x variables are (2, 11)
The standard deviations used for the beta priors are  (1.0444803217312628, 1.586983908089902)
The standard deviations used for the M12 priors are 0.09499285587637817

pm.model_to_graphviz(HAsfc81Model)

Check prior choice

with HAsfc81Model as model:
    prior_pred_HAsfc81 = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc81,dataZ["HAsfc81_z"].values,'HAsfc81')

Prior choice is as intended: Broad over the data range.

Sampling

with HAsfc81Model as model:
    trace_HAsfc81 = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [HAsfc81_M12_dist, HAsfc81_mu_M12, HAsfc81_b2_dist, HAsfc81_b2_beta, HAsfc81_b1_dist, HAsfc81_b0_dist, HAsfc81_sigmaY, HAsfc81_nuY, HAsfc81_sigma12, HAsfc81_mu_b2, HAsfc81_mu_b1, HAsfc81_mu_b0, HAsfc81_sigma2, HAsfc81_sigma1, HAsfc81_sigma0]

100.00% [15000/15000 11:57<00:00 Sampling 10 chains, 0 divergences]

Sampling 10 chains for 1_000 tune and 500 draw iterations (10_000 + 5_000 draws total) took 718 seconds.

with HAsfc81Model as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('HAsfc81'), 'wb') as buff:
            pickle.dump({'model': HAsfc81Model, 'trace': trace_HAsfc81}, buff)

Check sampling

with HAsfc81Model as model:
    dataTrace_HAsfc81 = az.from_pymc3(trace=trace_HAsfc81)

pm.summary(dataTrace_HAsfc81,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
HAsfc81_mu_b0	-0.06	0.81	-1.66	1.48	0.01	0.01	5438.0	3055.0	5456.0	4451.0	1.0
HAsfc81_mu_b1[0]	0.01	0.72	-1.40	1.39	0.01	0.01	4786.0	2800.0	4784.0	3776.0	1.0
HAsfc81_mu_b1[1]	-0.08	0.74	-1.54	1.35	0.01	0.01	4760.0	2996.0	4760.0	3738.0	1.0
HAsfc81_mu_b2[0]	-0.36	0.67	-1.69	0.90	0.01	0.01	3853.0	3574.0	3866.0	3483.0	1.0
HAsfc81_mu_b2[1]	-0.40	0.67	-1.68	0.92	0.01	0.01	4465.0	3899.0	4474.0	3963.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
HAsfc81_mu[275]	-0.16	0.17	-0.45	0.16	0.00	0.00	5195.0	5195.0	5325.0	4756.0	1.0
HAsfc81_mu[276]	-0.10	0.21	-0.44	0.34	0.00	0.00	4753.0	4660.0	5116.0	4932.0	1.0
HAsfc81_mu[277]	-0.16	0.17	-0.45	0.16	0.00	0.00	5195.0	5195.0	5325.0	4756.0	1.0
HAsfc81_mu[278]	-0.10	0.21	-0.44	0.34	0.00	0.00	4753.0	4660.0	5116.0	4932.0	1.0
HAsfc81_mu[279]	-0.16	0.17	-0.45	0.16	0.00	0.00	5195.0	5195.0	5325.0	4756.0	1.0

384 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc81,dataTrace_HAsfc81,'HAsfc81')

/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  warnings.warn(

with HAsfc81Model as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc81,'HAsfc81')

with HAsfc81Model as model:
    plotting_lib.pm.energyplot(trace_HAsfc81)

Posterior predictive distribution

with HAsfc81Model as model:
    posterior_pred_HAsfc81 = pm.sample_posterior_predictive(trace_HAsfc81,samples=numPredSamples,random_seed=random_seed)

/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  warnings.warn(

100.00% [2000/2000 00:03<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc81,posterior_pred_HAsfc81,dataZ["HAsfc81_z"].values,'HAsfc81')

plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc81,'HAsfc81',x1,x2)

Compare prior and posterior for model parameters

with HAsfc81Model as model:
    pm_data_HAsfc81 = az.from_pymc3(trace=trace_HAsfc81,prior=prior_pred_HAsfc81,posterior_predictive=posterior_pred_HAsfc81)

arviz.data.io_pymc3 - WARNING - posterior predictive variable HAsfc81_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc81,'HAsfc81')

plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc81,'HAsfc81',x1,x2)

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc81,'HAsfc81')

df_hdi_HAsfc81 = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc81,'HAsfc81',x1,x2)

df_hdi_HAsfc81

	Treatment_i	Treatment_j	hdi_ConfoMap_2.5%	hdi_ConfoMap_97.5%	isSignificant_on_ConfoMap	hdi_Toothfrax_2.5%	hdi_Toothfrax_97.5%	isSignificant_on_Toothfrax
0	Dry grass	Dry bamboo	-0.009401	0.101134	False	0.018900	0.136196	True
1	Dry lucerne	Dry bamboo	-0.030370	0.087266	False	-0.001664	0.109981	False
2	Dry lucerne	Dry grass	-0.073950	0.050741	False	-0.087209	0.039690	False
3	BrushNoDirt	BrushDirt	-0.518195	0.206099	False	-0.147466	0.286371	False
4	Control	BrushDirt	-0.376578	0.303712	False	-0.030208	0.268768	False
5	Control	BrushNoDirt	-0.097626	0.347835	False	-0.157188	0.295603	False
6	RubDirt	BrushDirt	-0.493950	0.236296	False	-0.128138	0.198331	False
7	RubDirt	BrushNoDirt	-0.231775	0.307111	False	-0.273242	0.218386	False
8	RubDirt	Control	-0.321828	0.132966	False	-0.274491	0.088330	False
9	Clover+dust	Clover	-0.455539	0.323286	False	-0.516906	0.156140	False
10	Grass	Clover	-0.656279	0.348255	False	-0.691228	0.029980	False
11	Grass	Clover+dust	-0.454243	0.139046	False	-0.317566	0.103036	False
12	Grass+dust	Clover	-1.007356	-0.112084	True	-0.965316	-0.269417	True
13	Grass+dust	Clover+dust	-0.689484	-0.292849	True	-0.567226	-0.229459	True
14	Grass+dust	Grass	-0.708570	-0.038109	True	-0.516849	-0.061466	True

if writeOut:
    df_hdi_HAsfc81.to_csv(outPathData+ 'hdi_{}.csv'.format('HAsfc81'))

Bimodal distribution in contrast plots of HAsfc81

For e.g. the pair Clover+Dust vs. Clover shows an unexpected bimodal distribution of the contrast. We now examine the traces carefully to exclude errors in the sampling:
Get the traces on ConfoMap of the interactions

m12_clover = trace_HAsfc81['HAsfc81_M12'][:,0,2]
m12_cloverDust = trace_HAsfc81['HAsfc81_M12'][:,0,3]
diff_m12 = m12_cloverDust-m12_clover

Get the traces of the treatments:

b2_clover = trace_HAsfc81['HAsfc81_b2'][:,2]
b2_cloverDust = trace_HAsfc81['HAsfc81_b2'][:,3]
diff_b2 = b2_cloverDust-b2_clover

Look at all the pairs

sns.pairplot(data=pd.DataFrame.from_dict(
    {'m12_clover':m12_clover,'m12_cloverDust':m12_cloverDust,'diff_m12':diff_m12,
     'b2_clover':b2_clover,'b2_cloverDust':b2_cloverDust,'diff_b2':diff_b2
    }
));

Let's have that of the differences again

sns.jointplot(x=diff_b2,y=diff_m12);

We see two sets that are distributed along parallel lines in the scatter plot. This means that the model estimates two subsets of possible differences.
However, when looking at the raw data at 'R_analysis/plots/SSFA_Sheeps_plot.pdf' one can see that the distribution of values for HAsfc81 on Clover (Sheep) measured by ConfoMap appear to have a bimodal distribution.
Thus, in combination with the chosen uninformative priors, the model correctly describes the distribution as bimodal.
In summary, we see no issue with the modeling and sampling.

Summary

Set the surface parameters for every treatment dataframe:

df_hdi_Asfc["SurfaceParameter"] = "Asfc"
df_hdi_HAsfc9["SurfaceParameter"] = "HAsfc9"
df_hdi_HAsfc81["SurfaceParameter"] = "HAsfc81"
df_hdi_R["SurfaceParameter"] = "R²"
df_hdi_epLsar["SurfaceParameter"] = "epLsar"

df_hdi_total = pd.concat([df_hdi_epLsar,df_hdi_R,df_hdi_Asfc,df_hdi_HAsfc9,df_hdi_HAsfc81],ignore_index=True)

Show the treatment pairs and surface parameters where the softwares differ

df_summary = df_hdi_total[df_hdi_total.isSignificant_on_ConfoMap != df_hdi_total.isSignificant_on_Toothfrax][["Treatment_i","Treatment_j","SurfaceParameter","isSignificant_on_ConfoMap","isSignificant_on_Toothfrax","hdi_ConfoMap_2.5%","hdi_ConfoMap_97.5%","hdi_Toothfrax_2.5%","hdi_Toothfrax_97.5%"]]
df_summary

	Treatment_i	Treatment_j	SurfaceParameter	isSignificant_on_ConfoMap	isSignificant_on_Toothfrax	hdi_ConfoMap_2.5%	hdi_ConfoMap_97.5%	hdi_Toothfrax_2.5%	hdi_Toothfrax_97.5%
2	Dry lucerne	Dry grass	epLsar	True	False	0.000228	0.001524	-0.000167	0.001116
10	Grass	Clover	epLsar	False	True	-0.000153	0.002973	0.000293	0.003795
18	BrushNoDirt	BrushDirt	R²	True	False	-0.002096	-0.000672	-0.000420	0.000635
22	RubDirt	BrushNoDirt	R²	True	False	0.000181	0.002386	-0.000556	0.000619
24	Clover+dust	Clover	R²	True	False	-0.001005	-0.000290	-0.000138	0.000401
26	Grass	Clover+dust	R²	False	True	-0.000639	0.000580	-0.000792	-0.000199
28	Grass+dust	Clover+dust	R²	False	True	-0.001314	0.000824	-0.000515	-0.000045
33	BrushNoDirt	BrushDirt	Asfc	True	False	0.013059	7.915624	-2.631889	6.887819
35	Control	BrushNoDirt	Asfc	True	False	0.022112	14.415237	-0.616464	13.182578
36	RubDirt	BrushDirt	Asfc	False	True	-0.259987	20.831039	1.502426	13.522475
60	Dry grass	Dry bamboo	HAsfc81	False	True	-0.009401	0.101134	0.018900	0.136196

if writeOut:
    df_summary.to_csv(outPathData+ 'summary.csv')

Write out

!jupyter nbconvert --to html Statistical_Model_TwoFactor.ipynb

[NbConvertApp] Converting notebook Statistical_Model_TwoFactor.ipynb to html
[NbConvertApp] Writing 24528302 bytes to Statistical_Model_TwoFactor.html

!jupyter nbconvert --to markdown Statistical_Model_TwoFactor.ipynb

[NbConvertApp] Converting notebook Statistical_Model_TwoFactor.ipynb to markdown
[NbConvertApp] Support files will be in Statistical_Model_TwoFactor_files/
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Making directory Statistical_Model_TwoFactor_files
[NbConvertApp] Writing 109106 bytes to Statistical_Model_TwoFactor.md